#include<cstdio>
#include<iostream>
#include"stick.h"
using namespace std;

int main()
{
	// test for sphere
	Parameters par;
	par.for1.a = 1; // parameters for function surface1 which describe first
	par.for1.b = 1; // stick surface. Surface should be differentiable in any
	par.for1.c = 1; // points; Of course numbers of parameters for first stick
	par.for1.d = -1; // can be different than there - only redefine Parameters

	par.for2.a = 1; // parameters for function surface2 which describe second
	par.for2.b = 1; // stick surface. Position of center of first stick should
	par.for2.c = 1; // be always (0,0,0), but start point for center of second
	par.for2.d = -1; // stick should be different than (0,0,0);
	
	par.forPhi.precision = 100; // how many steps we want to do for calcule
   	// area of intersection of excluded space;
	par.forPhi.current = 0; // first step of calculation

	par.zPlain = 1.5; // value of x axe which define where a plain which
		// cut the intersection is.
	
	const int n=6; // number of variables
	double x0[n] = {0.8,0.,0.75, 1.4, 0, 1.5}; // start point. First 3 values are
		// aproximate position of common point for surface of first and second
		// sticks. Last three points are aproximate position of center of
		// second stick. Look that variable z should be equal to par.zPlain
		// and variable x should be equal to cos(2*PI*current/pricision) times
		// aproximate radius of intersection. Similar for y but with sin.

	double R = 1;
	double *points[n];
	for(int i=0; i<n; i++) points[i] = new double[par.forPhi.precision];
	
	useGSL(&function, (void*)&par, n, x0, 0, points);
	for(int i=1; i<par.forPhi.precision; i++)	
	{
		par.forPhi.current = i;
		for(int j=0; j<n; j++) { x0[j] = points[j][i-1]; }
		x0[0] = R*cos(2*M_PI*par.forPhi.current/par.forPhi.precision)+1E-10;
		x0[1] = R*sin(2*M_PI*par.forPhi.current/par.forPhi.precision)+1E-10;
		x0[3] = 2*R*cos(2*M_PI*par.forPhi.current/par.forPhi.precision)+1E-10;
		x0[4] = 2*R*sin(2*M_PI*par.forPhi.current/par.forPhi.precision)+1E-10;
		useGSL(&function, (void*)&par, n, x0, i, points);
	}

	printf("\nIntersection area is amroximal equal to %.6e\n",\
		calculateArea(points[3], points[4], par.forPhi.precision));

	for(int i=0; i<2; i++) delete[] points[i];
	return 0;
}
